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\o ■ 

("■^ ■ The aim of this paper is to review and discuss qualitatively some results on the properties of 

(— ^ ' amorphous packings of hard spheres that were recently obtained by means of the replica method. 

rvQ I The theory gives predictions for the equation of state of the glass, the complexity of the metastable 

states, the scaling of the pressure close to jamming, the coordination of the packing and the pair 
correlation function in any space dimension d. The predictions compare very well with numerical 
simulations in d = 2, 3. The asymptotic predictions for d — > 00 are within the rigorous bounds. The 
theory can be extended to binary mixtures and to hard core potentials with an attractive tail or 
square well. 



X. 



a 



The study of amorphous packings of hard spheres is relevant for a large class of problems, including liquids, glasses, 

(— I ' colloidal dispersions, granular matter, powders, porous media, etc. and a large amount of numerical and experimental 

(~| ] data is available in the literature [1-17]. Moreover, the problem of sphere packing is related to many mathematical 

' . problems and arises in the context of signal coding and of error correcting codes, therefore it has been investigated in 

■ ^H ' detail by the information theory community [18, 19]. Nevertheless, a satisfactory characterization of the amorphous 

^ , states of a system of identical hard spheres is not yet available and in particular the question whether a glass transition 

exists in finite dimension is still open. From a rigorous point of view, for space dimension d > 3 only some not very 

restrictive bounds on the maximal density have been obtained, and in particular it is still unclear whether the densest 

H ■ packings for d — > 00 are amorphous or crystalline. 

i-rt [ Recently, a quantitative description of the glass transition in structural glasses has been obtained by means of the 

^ i replica trick [20-25]. This method was successfully applied to Lennard- Jones systems [20-23] and, more recently, to 

O ' hard spheres in R [24-26], on the Bethe lattice [27], and on the hypercube [28]. In this approximation the glass 

1^1 transition turns out to be similar to the 1-step replica symmetry breaking (IRSB) transition that happens in a class 

of mean-field spin glass models and indeed the replica strategy described above was inspired by the exact solution of 

CN . these models [29-32]. 

Jft I In this paper I review these results, for the case of hard spheres in M , without discussing the replica method, for 

L^ which the reader can refer to the original papers [20-22, 25, 26]: here I will discuss in details the limits of the theory 

\^ and how it compares with numerical results. As a matter of fact, despite the strong idealizations involved in the 

t:;^ 1 theory, the agreement with numerical data is surprisingly good. 

^: 

O . I. WHAT IS AN "AMORPHOUS PACKING"? 

a- . . . 

a In this section we will try to define "amorphous packings" of hard sphere systems and in particular the "random 

', close packing" density, i.e. the maximum density of amorphous packings. We will see that there are difficulties 

i-rt in the theoretical definition of this concept that allow only to define it approximately. The usual way to estimate 

P5 ^ this quantity in experiments or numerical simulation is to compress the system according to some protocol: e.g. in 

O ■ numerical simulations the particle diameter is slowly increased during a molecular dynamics run, or in experiments 

^ ' particle are thrown randomly in a box and the box is subsequently shaked [1-15]. The density of the system increases 

^ , with time and usually approaches a value of fRCP ~ 0.64. However the precise value of the latter quantity is found 

l^ '■ to depend in a non trivial way on the details of the experimental protocol. Then, to find a theoretical estimate of 

'^jj ! V-RCP we will have to refer to some idealization. For concreteness we will discuss the problem in M but the discussion 

j^ applies to any space dimension. 
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FIG. 1: Schematic phase diagram of hard spheres in M . (Left) Continuation of the liquid equation of state in the metastable 
region. (Right) Expected behavior in presence of a thermodynamic glass transition. 



A. The fate of the liquid above freezing: a naive definition of random close packing 



The equilibrium phase diagram of hard spheres in R is sketched in Fig. 1, where the density is reported as a function 
of the 'packing fraction ip = pVd{D/2), D being the diameter of a sphere and Vd{r) the volume of a d-dimensional 
sphere of radius r, so that ip is the fraction of volume covered by the spheres. At low density the system is in a liquid 
phase (as defined e.g. by the low density virial expansion); the maximum possible density is realized by the FCC 
lattice. A first order phase transition between the liquid phase and the FCC crystal phase is found by numerical 
simulations. 

The first naive idea to define amorphous states of hard spheres at high density is to assume that the liquid phase 
can be continued above the freezing density ipf and to look at its properties at large density. For instance one can 
choose a functional form that represents well the equation of state of the liquid below iff (e.g. the Carnahan-Starling 
or Percus-Yevick equation of state [33]) and assume that it describes the liquid phase also above ipf. On increasing 
the density the pressure of the liquid increases as the average distance between particles decreases: one may expect 
that it diverges at a point where the particles get in contact with their neighbors, and the system cannot be further 
compressed. Then one can call this point random close packing density (RCP), see left panel of Fig. 1. 

The first objection to this proposal is that an intrinsic stability limit of the liquid (a spinodal point) might exist at 
a density above </?/ due to thermodynamic or kinetic reasons. A thermodynamic spinodal is likely not to be there and 
any reasonable continuation of the liquid equation of state do not predict such an instability (manifested e.g. by an 
infinite compressibility). A kinetic spinodal, related to the existence of the crystal [34], could instead exist, at least in 
monodisperse systems. This would imply the impossibility of reaching amorphous jammed states if the compression 
rate is not very high. We will assume in the following that the metastable liquid can be compressed as slow as wished 
avoiding crystallization. This is not a very good assumption for monodisperse systems but seems to be more close to 
reality for a suitably chosen binary mixtures. This point requires a better investigation, e.g. following the analysis of 
[34], and we will not discuss it further here. 

A second objection is that in presence of a first order phase transition the continuation of one phase into its 
metastable region is not well defined due to the appearance of essential singularities at (pf . Many possible continuations 
of the low density equation of state above ipf are possible: this means that the properties of the "liquid" above ipf 
will depend on the history of the sample, much as it happens for the hysteresis of a ferromagnetic system. This 
means that, even neglecting crystallization, the definition of RCP above is not precise from a theoretical point of 
view. Nevertheless, the ambiguity is expected to be exponentially small in the distance from iff, and as the distance 
between iff and the maximum density fpcc is not so large, one might expect to obtain a meaningful result in this 
way. And indeed the difference between different possible continuation of the liquid equation of state (e.g. Carnahan- 
Starling, Percus-Yevick, Hypernetted Chain) is of the order of 10% above (pf. 

The main problem with this definition is that all these equations predict a divergence of the pressure at unphysical 
large values of (p: e.g. the CS equation predicts (pRCP = 1 which is clearly wrong since it is larger than the FCC value 
and moreover implies that the available volume is completely covered by the spheres. Thus this idea is definitely too 
naive and we have to resort to something more refined. 



B. The ideal glass transition: a (still naive) definition of RCP 

A way out of this contradiction is to assume the existence of a thermodynamic glass transition in the metastable 
hquid branch of the phase diagram. Such a transition can be expected for various reasons but we will not try to 
justify it: we will simply assume that it is there and investigate the consequences of such an assumption. 

So let us assume again that the liquid phase can be continued above (/?/ and neglect the (small) ambiguity in 
its definition due to its intrinsic metastability with respect to crystallization. We assume that at a density ipn a 
thermodynamic glass transition (sometimes called ideal glass transition) happens. The transition is signaled by a 
jump in the compressibility of the system. A simple qualitative argument to explain this is the following: in the 
dense liquid phase particles vibrate on a fast time scale in the cages made by their neighbors, while on a much larger 
time large scale cooperative relaxation processes happen (structural relaxation). If we change the density by Aip the 
pressure will instantaneously increase by a APq; very rapidly the average size of the cages will decrease a little due 
to the increase in density and the pressure will relax to a value APy < APq. Then, on the time scale of structural 
relaxation, the structure will change to follow the change in density and the pressure will relax further to a value 
APoo < AP^. At the glass transition the latter relaxation is frozen as the corresponding time scale becomes infinite. 
Thus, in the glass phase the increase in pressure following a change in density will be larger than in the liquid phase, 
leading to a smaller compressibility K = (p^^{Aip/AP). 

The schematic phase diagram that we expect in presence of a glass transition is in right panel of Fig. 1. It is evident 
that the existence of a glass transition can cure the paradoxical behavior of the pressure of the liquid, that seems 
to diverge at a density bigger than the FCC density. Then we can define the RCP density as the point where the 
pressure of the glass diverges. 

Still this definition of RCP in unsatisfactory for two reasons. The first is that the ambiguity in the definition of the 
liquid equation of state due to its metastability affects also the glass: so the glass equation of state will be theoretically 
not well defined, with an ambiguity of the order of 10%, depending on the equation of state one chooses to describe the 
liquid. Moreover, the existence of the crystal means that one can construct configurations of the system representing 
mixtures of the crystal and the glass states; such configurations will have an arbitrary degree of local order and will 
span the whole range of densities between the density of the glass and that of the crystal: explicit examples have 
been constructed in [12]. It is very difficult to define an "order metric" to quantify the order in a given configuration 
of a large system. Thus on a practical ground it is very difficult to distinguish a "pure glass state" from a mixture of 
glass and crystal. To resume, the existence of the crystal poses both theoretical and practical difficulties in defining 
RCP as the point where the pressure of the equilibrium glass diverges. 

A second difficulty is the following. The glass transition, in the standard picture coming from the analysis of mean 
field models, is due to the appearance of many metastable states: in the case of hard spheres this means that we 
expect that, in addition to the stable glass states corresponding to the equilibrium glass, we will have many metastable 
packings. These packings can be described as follows: in the equilibrium liquid at a given (high) density ip < i^k 
the particles will vibrate around local stable structures, that are visited subsequently on the scale of the structural 
relaxation. If we now "artificially" froze the structural relaxation and compress the system, the pressure will increase 
faster than if the structure is allowed to relax for the same reason as above: the system will be forced to reduce the 
size of the cages to respond to a change in density. Then the pressure will diverge at the point where the particles 
get in contact with their neighbors and the average size of the cages is zero. In this way we can produce jammed 
configurations in a whole range of densities tp < ^rcp [12-15]. Experimentally or numerically this can be done by 
compressing fast enough in order to disallow the system to relax the structure during compression. And it is clear 
that, as the time scale for structural relaxation is expected to diverge on approaching ipx, at some point it will fall 
beyond any experimentally accessible time scale and necessarily the system will be frozen into a metastable state. 
Thus the ideal glass states are likely to be unobservable in practice: one will observe instead a nonequilibrium glass 
transition to a state that, again, will depend on the experimental protocol. 

To summarize, the difficulties in the definition of the random close packing arc mainly due to metastability, either 
of metastable glass states with respect to the ideal glass states, or of the whole glass states with respect to the crystal. 
Thus, in the following we will assume the idealized phase diagram given in Fig. 1, neglect the existence of the crystal 
and try to compute properties of the pure ideal glass state, and possibly also of metastable glass states. We will see 
that even in this "zeroth order" approximation we can obtain a very good agreement with numerical simulations. 
Moreover, we will consider the limit d — > oo in which the ambiguities related to metastability are expected to be much 
smaller. We believe that it is possible that in this limit the theory is exact. 

The method we used has been described in [25, 26]. In an extended version of this paper we will present a simpler 
derivation, explain how crystallization is removed in the theory and qualitatively discuss the effect of metastability. 
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FIG. 2: (Left) Equation of state in d = 3 using the Carnahan-Staling equation of state to describe the liquid. The inverse of 
the reduced pressure PP/ p is reported as a function of ip. The ideal glass transition density is ipK = 0.62 and the random close 
packing density is ipRCP = 0.68. The metastable state reported in the figure correspond to a jamming density tpj = 0.658. 
Numerical data from [10] are reported for comparison. (Inset) Complexity as a function of the jamming density for d = 3. 
(Right) Scaled correlation function gj(r)/gj{D) close to r = D. The theory is compared with numerical data from [14]. 

C. On the protocol dependence of the random close packing density 

Generally speaking, we can use a given algorithm to produce jammed configtirations of hard spheres: the case we 
discussed above corresponds, for instance, to a molecular dynamics simulation during which the system is compressed 
at a given rate. The algorithm will stop when the system is jammed: the final density will be a random variable 
depending on the initial data and possibly on some randomness built in the algorithm itself. We can define the 
probability V{ifj) of reaching a final density ipj: it will be the product of the density of states p{ipj), i.e. the 
number of jammed configurations with density (pj, times the probability /3(</?j) that the algorithm finds a particular 
configuration with density ip>j. The former quantity depends only on geometrical properties of the configuration space 
of hard spheres, while the latter encodes the properties of the algorithm [15]. 

We can compute the complexity T,{ipj) = hmAr^oo J^^^^og p{ipj) for the "amorphous states" defined above. This 
quantity is reported in Fig. 2. It vanishes linearly at fRCP giving an exponential distribution of the jammed states, 
Pifj) "^ exp[A^(</5j — (pRCp)]- However the quantity f3{ipj) is also expected to scale as P{<pj) ~ exp[Nb{(pj)]; the 
resulting V{ipj) will be strongly peaked around a given value ipj - the maximum of S((/3j) + h{(fj) - [35], which 
however might depend on the particular algorithm through the function f3{ipj), giving rise to the protocol dependence 
observed in experiments, see [15] for a detailed discussion of this important issue. However, the value fRCP defined 
as the point where 'S(ipj) vanishes is a property of the system that does not depend on the algorithm (at least if one 
accepts the idealizations discussed in the previous subsection, i.e. if one neglects the existence of the crystal). 



II. RESULTS 



A. Equation of state of the glass 



The only input of the theory is the entropy of the liquid, S{(p). As discussed above, this quantity is not well 
defined above the freezing density and we have to choose an expression which describes well the liquid below (pf and 
extrapolate above ipf. 

d^l - In dimension rf = 1 we do not expect any glass transition. The entropy can be computed exactly and is given 
by S{p) = 1 — logp + log(l — pD): it diverges at the close packing pD ~ 1. Indeed we find no glass transition. 

d^2 - In d = 2 we can use either the Henderson expression for S{ip) [36], or the improved expression of Luding [37], 
with very small quantitative differences. Again we find no glass transition, consistently with the numerical simulations 
of [38], but we find a glassy solution which never becomes stable. The properties of the amorphous states seem to be 
similar to the one observed in the numerical simulations of [39] (see below) . However in d = 2 the hard sphere liquid 
has very peculiar properties so this case deserves further investigation. 



d^3 - In d = 3 we used the Carnahan-Starling expression [33] for the entropy, which reproduces very well the 
numerical data for the equation of state. We predict a glass transition density p^ ^ 0.62 in good agreement with 
numerical results. The pressure of the equilibrium glass diverges at fRCP ^ 0.68. We can also compute the equation 
of state of the metastable (or nonequilibrium) glasses. The latter are labeled by their jamming density ipj which is 
the density at which the pressure of the metastable glass diverges. We find a whole distribution of metastable glasses 
with jamming densities ipj < (pRCP- A comparison with numerical data is in Fig. 2: the comparison suggests that the 
value (fj ~ 0.64 observed in numerical simulation correspond to the jamming density of a metastable state in which 
the system is frozen due to the finite accessible time scale. The ideal value of (pRCP ~ 0.68 would be reached only by 
an infinitely slow compression (if crystallization is avoided) . 

d^ oo - The entropy of the hard sphere liquid for d ^ oo was computed in [40, 41] in two different ways, and in both 
cases it was found that S{(p) is given by the ideal gas term plus the first virial correction (i.e. by the Van der Waals 
equation), S{ip) = 1 — log/?— 2'^^^(p with exponentially small correction in d ^ oo. Using this expression we find that 
both (fiK and (fiRCP scale asymptotically as d\ogd/2'^, with (ppcp — 'Pk ~ d/2'^. We can compare this prediction for 
the maximum density of amorphous packings with the best available bounds on the density of crystalline packings. 
Unfortunately, the best lower bound for periodic packings is the Ball bound ip > d/2'^, while the best upper bound is 
the Blichfeldt's one, ip < 2~'^'^ [18, 19]. Our result for prcp lies between these bounds so we cannot give an answer 
to the question whether the densest packings of hard spheres in large d are amorphous or crystalline. Hopefully better 
bounds on the density of crystalline packings will address this question in the future. Recent related work appeared 
also in [42]. 

B. Scaling close to jamming 

We can derive asymptotic relations for the behavior of the metastable states at </? ^ </?j (the equilibrium glass 
corresponds to the particular case pj = 'Prcp)- 

The complexity S((y9j), that represents the logarithm (per particle) of the number of jammed structures with density 
(/?j, is given by 



nv.j)^s{^.j)-d\og 



V8 



2dpjY{ipj) 



f, (1) 



where Y{pj) is the value of the pair correlation at contact. It is reported in Fig. 2 for d = 3. Note that unfortunately 
we cannot investigate the stability of the metastable packings so we do not know the minimal value of pj. 
The entropy of a metastable glass is found to diverge as 

s{p, pj) - d\og{pj - p) (2) 

close to Pj, so that the (reduced) pressure diverges as 

ds{p,pj) dpj 

p[p, pj) = -p . (3) 

dp Pj - p 

It is important to remark that the corrections to the leading order seem to be quite large, as observed in numerical 
simulations. 

C. Correlation function 

The most interesting results on the scaling close to pj are obtained by investigating the pair correlation function 
gj{r) of the spheres in the packing. This quantity can be (partially) computed within our theory: the details are in 
[25]. We find that gj{r) = g{r)[l + Sgj{r)], where g(r) is the pair correlation of the liquid and Sgj{r) is a correction 
which is large for r — D < ^p — pj. In this region we have, for p ^* pj: 

/tF p((^, pj){r — D 






2 D 



Ao{u) = 2 dyye-y -"^""y = l-VT^ue'' {l-eTi{u)) . (4) 



In Fig. 2 we report numerical data from [14] on the delta-peak contribution for a packing in d = 3 with pj ~ 0.64 
and Pj — p ^ 10""'^^. We see that our result, given by Eq. (4) for pj = 0.64, perfectly agrees with the data. 
Unfortunately there are some contributions to gj{r) that we still are not able to compute: we expect these corrections 
to be responsible for the power law divergence of gj (r) for r ^ D and for the splitting of the second peak as described 
e.g. in [14]. 
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FIG. 3: Number of contacts. (Left) Data from [14] in d = 3 and </pj ~ 0.64. The theory predicts z = 6 + 20y'(pj — ip. (Right) 
Data from [39] in d = 2 and (fj ~ 0.83. The theory gives 2 = 4+ 14.5^/ipj — tp. 



D. Number of contacts 



For ip — > ifij there is a shell of width y^(p,j — f around a given particle where the probability of finding other 
particles is very high, as expressed by the correction Sgj{r). These particles will become neighbors of the particle in 
the origin at ip = ipj, and their number is finite as the integral of 5gj{r) on this shell is finite. The integral of gj{r) 
on a shell D < r < D + 0{y^(pj — ip) gives then the average number of contacts for (p — > ipj, which turns out to 
be z = 2d[l + A((pj)yyipj — ip\, where the coefficient A(ipj) can be explicitly computed. For ifj — ip finite but small 
this number can be interpreted as the number of particles that collide with the particle in the origin during a finite 
but long time t [14, 39]. The result is compared with numerical data in Fig. 3. Remarkably it compares well with 
simulations also in d = 2, indicating that the (unstable) glass phase we find in d := 2 might have some connection 
with bidimensional amorphous packings. 



III. CONCLUSIONS 

The replica method seems to be a powerful tool to investigate the properties of amorphous packings of hard spheres. 
Despite the idealizations involved in this approach, the results seem to compare well with numerical simulations in 
d = 2,3. As the method works in any space dimension, it gives predictions in d > 3 that we hope will be tested in 
the future. In the limit d — > oo the method is likely to become exact and we hope that this work will stimulate new 
more rigorous investigations of this (and related) problems. 

There are still a number of issues that have to be addressed within the theory: for instance, a full computation 
of the pair correlation of the glass is mandatory, and an investigation of the stability of the jammed configurations 
should be performed. 

Work is in progress to extend the method to the case of binary mixtures, following [23]. In this case crystallization 
should be less relevant and more reliable prediction should be obtained, as in the case of Lennard- Jones systems. 
We would like also to apply the method to hard spheres with attractive potentials like Yukawa tails or square wells: 
for these potentials Mode-Coupling theory predicts a rich phenomenology (reentrance of the glass transition line, 
glass-to-glass transition, ...) [43]. Some of these results have indeed already been reproduced within a replica method 
similar to the one presented here [44]. 
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